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ABSTRACT 

Recent gamma-ray observations of middle-aged supernova remnants revealed a mysterious broken power-law 
spectrum. Using three-dimensional magnetohydrodynamics simulations, we show that the interaction between 
a supernova blast wave and interstellar clouds formed by thermal instability generates multiple reflected shocks. 
The typical Mach numbers of the reflected shocks are shown to be M 2 depending on the density contrast be- 
tween the diffuse intercloud gas and clouds. These secondary shocks can further energize cosmic-ray particles 
originally accelerated at the blast-wave shock. This "two-step" acceleration scenario reproduces the observed 
gamma-ray spectrum and predicts the high-energy spectral index ranging approximately from 3 to 4. 

Subject headings: ISM: supernova remnants — magnetic fields — turbulence — acceleration of particles 



1. INTRODUCTION 

Supernova remnants (SNRs) are believed to be the sites of 
Galactic cosmic-ray acceleration through a diffusive shock 
acceleration mechanism (DSA; see, e.g., Blandford & Eich- 
ler 1987) and multi- wavelength nonthermal emissions from 
SNRs originating in accelerated particles have been detected 
(Koyama et al. 1995, Aharonian et al. 2008a). Recent obser- 
vations using Fermi Gamma-ray Space Telescope revealed the 
gamma-ray spectra from middle-aged SNRs; W51C, W28, 
and W44 (Abdo et al. 2009b; 2010c; 2010d). The observed 
gamma-ray spectra is likely to be explained by accelerated nu- 
clei whose energy distributions are universally described by 
a broken power-law with an break energy of approximately 
10 GeV These SNRs are also known to be interacting with 
molecular clouds. In a neutral medium such as in the vicinity 
of molecular clouds, the damping of magnetohydrodynamic 
(MHD) waves by ion-neutral collisions leads to escape of 
particles from the shock wave that interrupts acceleration at 
some break energy (Abdo et al. 2009b). Theoretical studies 
of the particle acceleration that takes into account the effect 
of the wave dumping well explain the observed break energy 
of approximately 10 GeV for middle-aged SNRs (Ptuskin & 
Zirakashvili 2003). In order to understand the observed bro- 
ken power-law spectrum, however, we should still find some 
mechanisms that accelerate particles beyond the break energy 
whose spectrum continues, at least, to the TeV energy level 
as a power-law. The power-law index s of the accelerated 
particles above the break point is roughly ^ ~ 3-4, which is 
much steeper than the conventional DSA scenario that pre- 
dicts s 0:^2. So far several scenarios have been proposed to 
explain the observed spectra (Malkov et al. 2010; Ohira et 
al. 2010b; Uchiyama et al. 2010; Li & Chen 2010). In this 
Letter, we propose another one. 

2. CLOUDY INTERSTELLAR MEDIUM 

To Study the physics of the shock-cloud interaction in 
SNRs, we examine the propagation of a shock wave through 
a medium generated as a natural consequence of the thermal 
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instability. The thermal instability is a mechanism for conden- 
sation of interstellar gas driven by runaway radiative cooling 
(Field 1965). Theoretically, the thermal instability is shown 
to be effective during the cloud formation process (Koyama 
& Inutsuka 2002; Hennebelle et al. 2008; Inoue & Inutsuka 
2008; 2009), and it always dominates the dynamics of the 
cloud formation (Heitsch et al. 2008). In this work, the three- 
dimensional, ideal magnetohydrodynamics equations includ- 
ing the effects of radiative cooling, heating, and thermal con- 
duction are solved using the Godunov-CMOC-CT scheme 
(Inoue et al. 2009; Sano et al. 1999; Clarke 1996). We 
perform simulations in (2 pc)^ numerical domain with a grid 
spacing Ax = 2/1024 pc 2 x 10~^ pc, imposing a periodic 
boundary condition. With this resolution, we can safely ex- 
press the smallest structure formed by the thermal instability 
(Koyama & Inutsuka 2004; Inoue et al. 2006). Panel (a) of 
Fig. [T] shows the result of the cloud formation by the thermal 
instability (detailed dynamics of the cloud formation can be 
found in Inoue & Inutsuka 2008; 2009; Inoue et al. 2009). 
The regions in blue represent the condensed clouds (ric ^ 40 
cm~^, Tc 90 K), while the diffuse intercloud gas is depicted 
in yellow (rii c^O.l cm~^, 7] 5000 K). We have initially im- 
posed the j-directional uniform magnetic field with strength 
B = 5 /iG, which is typical in the interstellar medium (Bech 
2000). Since the most unstable scale of the thermal instability 
is ~ 1 pc and the clouds are formed by condensation along 
the magnetic field line, the resulting clouds have sheet-like 
morphology whose scale is ~ 1 pc and thickness is ^0.1 
pc. Observationally, the structures formed by the thermal 
instability are expected in the envelope of molecular clouds 
(Sakamoto & Sunada 2003), where the interaction with a su- 
pernova shock wave takes place. 

3. FORMATION OF SNR 

By setting static, high-pressure plasma with P^/kB = 5.0 x 
10^ K cm~^ and = 0.1 cm~^ at the x = pc boundary, we 
inject a strong shock wave that propagates into the cloudy 
medium. The resulting average value of the shock velocity 
is 536 km s~\ which corresponds to that of the SNR with an 
age of ~ 10"^ yr. The density structure after the propagation 
of 3400 years is shown in the panel (b) of Fig.[T] In principle, 
a shock wave hitting a dense cloud should be stalled, because 
the momentum conservation results in a smaller velocity for 
medium with increased inertia. However, the shock wave 
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maintains high velocity in the diffuse intercloud medium that 
fills most of the volume. Thus, we can expect the conven- 
tional picture of DSA in most of the volume, even if the SNR 
interacts with dense clouds. 

Due to the velocity difference between the clouds and inter- 
cloud region, the shock is highly deformed. It is known that, 
even if the pre- shock medium is uniform, the curved shock 
wave leaves vorticity behind the shock wave (Crocco's theo- 
rem), driving turbulence in the SNR. As shown in the panel (c) 
of Fig.[TJ in such a turbulent downstream region, the magnetic 
field is amplified by the stretching of the field lines (Giacalone 
& Jokipii 2007; Inoue et al. 2009; see also Jones & Kang 1993 
who predicted the magnetic field amplification at the shear 
layer around the cloud caused by the shock-cloud interac- 
tion). In the present case, the average magnetic field strength 
in the SNR saturates approximately at ~ 40/iG, which is twice 
the maximum field strength achieved by the shock compres- 
sion alone. The maximum magnetic field strength saturates at 
~ 400/iG in the regions where local plasma (3 is on the order 
of unity. 

4. SECONDARY SHOCKS 

When a shock wave hits a dense object, a reflected shock 
wave forms and propagates back into the post-shock medium, 
in addition to the transmitted shock wave. In the present case, 
the pre-existing dense clouds in the pre- shock medium gener- 
ate a number of reflected shock waves. We call these reflected 
shock waves "secondary shocks" to contrast them with the 
original "primary shock" that propagates into the unshocked 
medium. Panel (d) of Fig. [T] shows the pressure distribution in 
the cross section of the SNR, indicating the existence of such 
secondary shock waves as discontinuous pressure jumps. The 
strength of the secondary shocks can be estimated by ana- 
lyzing the shock tube problem that describes the interaction 
between a dense gas and a pressure driven shock wave (Mi- 
esch & Zweibel 1994). The Mach number of the reflected 
shock wave Mr measured in the rest- frame of the shocked dif- 
fuse gas is obtained by solving a polynomial equation that is 
characterized by two parameters, namely the pressure ratio 

5 = Pc/Ps of the cloud and shocked gas and the density ratio 
e = ps/Pc of the shocked diffuse gas and the cloud (Miesch 

6 Zweibel 1994). Fortunately, the parameter S is very small 

10~^) in the present simulation. Expanding the polynomial 
with respect to (5, we find 
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ij+O(J) = 0, (1) 

where we have assumed the adiabatic index 7 = 5/3. Substi- 
tuting the typical density ratio e = 10~^ achieved in the simula- 
tion, we obtain Mr = 1.80. In the limit of e ^ (i.e., the very 
dense cloud case), we obtain Mr >/5. Because the primary 
shock is driven by isotropic thermal pressure and the transmit- 
ted shock wave is stalled heavily due to the large cloud den- 
sity, the primary shock tends to hit the cloud perpendicular 
to its surface. Thus, the above formula obtained by assuming 
one-dimensional geometry always enables us to evaluate Mr, 
even though clouds have complex structure. 

To show the strengths of the shock waves in the simulation, 
we plot the probability distribution function (PDF) of the pres- 
sure jumps in Fig.[2las a red line. Our algorithm to detect the 
pressure jump is as follows: We define the sphere whose ra- 
dius is 20 times numerical grid size and whose center is at 
each cell in the numerical domain. In each sphere, we calcu- 
late the minimum pressure Pmin, the maximum pressure Pmax, 
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Fig. 1. — Panel (a) shows a two-dimensional slice of density distribution 
of the cloudy medium as a consequence of thermal instability. Arrows in 
the panel indicate the orientation of the initial magnetic field. Panel (b) is 
a number density slice after 3,400 years shock propagation. Panel (c) is a 
slice of magnetic field strength at t = 3, 400 yr. Panel (d) is a pressure slice at 
t = 3, 400 yr. In the panels (a), (b), and (c), the x-y plane at z = pc are shown. 
In the panel (d), to display clearly the existence of the secondary shocks, the 
y-z plane at x = 1 .65 pc is shown. 

and the maximum pressure gradient, except for the injected 
hot plasma with ^ < 0.5 cm~^. If the maximum pressure gra- 
dient is larger than the critical pressure gradient Pmax/ (5 Ax), 
i.e., the pressure fluctuation in the sphere is caused within the 
narrow region, we regard a shock wave as being in the sphere 
whose pressure jump is AP = Pmax/Pmin- Here the factor 5 in 
the expression of the critical pressure gradient is chosen from 
the fact that the high-resolution shock-capturing scheme used 
in this study requires 5 grid points at most to express a shock 
wave. The bimodal distribution appearing in the PDF shows 
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Fig. 2. — Probability distribution function (PDF) of pressure jump AP. 
Red and blue lines are the results of the simulations with the pressure ratio 
6 ~ 10~^ and 10~^-^, respectively. Inner panel is the close up of AP G [1, 10]. 
Integrated area of the PDF in a specific range of the pressure jump is propor- 
tional to the total surface area of the shock waves whose strengths are in the 
given range. Thus, it is clear that the total surface areas of the secondary 
shocks is comparable to that of the primary shock. 

the existence of the primary and the secondary shocks. Ac- 
cording to the Rankine-Hugoniot relation, the Mach number 
of a shock for the gas of 7 = 5 /3 is related to the pressure 
jump as (Landau & Lifshitz 1959): 

M(AP)=^J^^^AP. (2) 

Substituting the peak of the PDF AP = 3.85 that corresponds 
to the secondary shocks, we obtain M = L81, which agrees 
well with the above analytic evaluation. To reinforce our anal- 
ysis, we have performed the additional simulation with larger 
primary shock velocity with Vsh = 2403 km s~\ i.e., the sim- 
ulation with smaller S (~ 10~^-^). Since Mr is almost inde- 
pendent of the parameter S, we should obtain almost the same 
PDF of the secondary shocks. The blue line in Fig. |2] shows 
the result of the additional simulation verifying our analytic 
estimate for the strengths of the secondary shocks. 

5. ORIGIN OF BROKEN POWER-LAW SPECTRUM OF COSMIC 
RAYS 

As mentioned earlier, recent gamma-ray observations of 
middle-aged SNRs have revealed a steep particle spectrum 
with power-law index s 3-4, while so:^2 below '^10 GeV. In 
the following, we consider DSA in the environment obtained 
by our simulation to explain the observed spectrum. Our idea 
is similar to that given in Jones & Kang (1993) who pointed 
out that the reflected shock wave caused by the shock-cloud 
interaction can affect the cosmic-ray acceleration in the SNR. 

The spectral index of the cosmic rays accelerated by the 
DSA mechanism at a shock wave with the Mach number M is 
given by (Blandford & Eichler 1987): 



Thus, the secondary shocks with M = 1.81 achieved in the 
simulation can accelerate particles with the spectral index 
^ = 3.75, while the primary shock with M ^ 10 gives ^ 2.0. 
In the galactic interstellar medium, the densities of the inter- 
cloud gas and HI clouds are determined by the balance of ra- 
diative cooling and heating. From the detailed evaluation of 
the heating/cooling rates, it is known that the density of inter- 
cloud gas is < 1 cm~^ and that of HI clouds is ^ 10 cm~^ 
between which gas is thermally unstable (see, e.g., Wolfire 



et al. 1995). Thus, the maximum value of the parameter e 
can be evaluated as emax = 4^i,max/^c,min — 0.4 that leads to 
the upper bound of the spectral index s ^ 4 from eqs. ([B, O 
and ([3]). If the preshock clouds are molecular and much denser 
than the HI cloud, the parameter e can be substantially smaller 
than 0.1 realized in the present simulation. From the typical 
Mach number of the secondary shocks in the dense cloud limit 
(M = >/5), the lower bound of the index would be 5- :^ 3. Thus, 
the spectral index of the particles accelerated at the secondary 
shock would be limited in the range 3 < ^ < 4, which agrees 
well with the recent observations of middle-aged SNRs (Abdo 
et al. 2009b; 2010c; 2010d). Note that recent numerical simu- 
lations have shown that the two-phase structure of clouds and 
diffuse intercloud gas is also generated even inside molecular 
clouds as a result of the thermal instability (see, e.g., Hen- 
nebelle et al. 2008). This suggests that our scenario would be 
unchanged even if supernova explosion happens inside molec- 
ular clouds. 

In contrast to the primary shock, the secondary shocks are 
located in the fully-ionized, postshock medium of the pri- 
mary shock. Thus, the acceleration process at the secondary 
shocks is free from the wave damping process due to the ion- 
neutral collisions, so that the acceleration beyond the break 
energy is possible. The maximum energy of accelerated pro- 
tons achieved at the secondary shock can be estimated using 
the DSA theory (Malkov & Drury 2001): 




where we have scaled using the average magnetic field 
strength in the simulated SNR ({\B\) ^40 /iG), and the degree 
of the magnetic field fluctuations r] = 6B^/B^ has evaluated us- 
ing the power spectrum of the magnetic field for the particles 
with energy near £^max. For the shock velocity, we have sub- 
stituted the primary shock velocity, since the sound speed in 
SNR, which is approximately the velocity of the secondary 
shocks, is always comparable to the primary shock speed 
(more precisely, Vsh,pri = 3cs/a/5 for the strong shock with 
adiabatic index 7 = 5/3). As for the acceleration timescale, 
we have substituted the age of SNR. The minimum lifetime 
of the secondary shocks is given by the radial sound cross- 
ing time of the SNR shell, which is typically one tenth of 
the age of SNR. Thus the maximum attainable energy esti- 
mated in eq. d?]) can be smaller by an order of magnitude, 
although the secondary shocks do not always propagate radi- 
ally. Note that the spatial extension of each secondary shock 
is roughly given by the length of the clouds, which is essen- 
tially determined by the most unstable scale of the thermal 
instability 1 pc. Thus, the scale of each secondary shock 
is much larger than the gyro-radius of the maximum energy 
proton /g = 0.017(£/10^'^eV)(5/40/iG)"^ pc, indicating that 
it does not restrict the maximum attainable energy. There- 
fore, in principle, the secondary shocks can accelerate parti- 
cles as large as 100 TeV in the middle-aged SNRs. Even if 
^max is an order of magnitude smaller owing to the lifetime 
of the secondary shocks, the maximum energy of 10 TeV is 
still large enough to account for the TeV gamma-ray emission 
from middle-aged SNRs (Abdo et al. 2009a; 2009b; 2010c; 
2010d; Aharonian et al. 2002; 2008b; Buckley et al. 1998). 

Since the Mach numbers of the secondary shocks are small, 
the injection rate of particles into the acceleration process in 
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the secondary shocks would be much smaller than that in the 
primary shock. However, suprathermal particles generated 
by the primary shock can be advected downstream and re- 
accelerated to higher energies at the secondary shocks. As 
shown in Fig.[2l the total surface area of the secondary shocks 
is comparable to that of the primary shock, suggesting that 
most advected particles meet the secondary shocks. If the vol- 
ume filling factor of the clouds is much larger than that real- 
ized in our simulations, the total surface area of the secondary 
shocks would be much larger. In that case the suprather- 
mal particles can be successively accelerated by multiple sec- 
ondary shocks. Given the particle momentum spectrum ac- 
celerated at the primary shock as No(p) oc exp(-/7//7br), 
where Phr — 10 GeV is the break momentum due to the wave 
damping, the spectrum re- accelerated n times by the sec- 
ondary shocks becomes (Melrose & Pope 1993): 

rp/d 

Nnip) = (s-l)(p/dr / q'-'Nn-i(q/d)dq, (5) 
Jo 



-2 

= P n^n 



OC 



Phvd 



2 (for p <C ;?br), 

' (for p > /7br), 



-]d^--^^^'\ (6) 
a J 

(7) 



where 



{gl)k--{gm)k r 
{hi\...{hl)k kV 



(8) 



is the generalized hypergeometric function written using the 
Pochhammer symbol [(g)k = g(g+l)...(g-^k-l), (g)o = 1], 
a = s-2, b = s-1, d = {(^- l)/(^ + 2)}^/^. Here we have as- 
sumed that the injection at the secondary shocks is inefficient 
owing to their small Mach number. Note that eq. ©is valid 
for pinj < p < Ejnsix/c, where pinj is the injection momentum 
above which particles can cross the shock. Moreover, we have 
assumed all the secondary shocks have the same Mach num- 
ber and thus the same acceleration slope s given by eq. ©. 
Thus, we can expect the unique broken power-law spectrum 
irrespective of the experienced number of the re- acceleration. 



6. PREDICTIONS AND DISCUSSIONS 

The above results allow us to make the following predic- 
tions, (a) The spectral index s above the break energy cannot 
be substantially smaller than 3, since the typical Mach num- 
ber of the reflected shocks in the dense cloud limit is Mr = >/5. 
Note that this is not true for young SNRs, because thermal 
pressure in the Sedov-phase shell decreases toward its cen- 
ter, which makes eq. ([B inapplicable. In that case, the Mach 
number of the secondary shocks would increase as they prop- 
agate toward the center of the SNR, leading to a shallower 
slope. The broken power-law spectra observed in a young 
SNR: Cassiopeia A (Abdo et al. 2010a) and possibly young 
SNR: IC443 (Abdo et al. 2010b) may be the case, (b) Mas- 
sive stars, progenitors of supernovae, are born in giant molec- 
ular clouds. Thus, most of middle-aged SNRs in the galactic 
mid-plane will interact with molecular clouds, and will have 
a broken power-law cosmic-ray spectrum whose index above 
the break energy is bounded by 3 < ^ < 4. (c) If the spectral 
index above the break energy at a middle-aged SNR is sub- 
stantially smaller than 3, the softening of the spectrum would 
be attributed to particle escape. In that case the spectral index 
substantially depends on the evolution of maximum attain- 
able energy £^max(0 and injection rate r](t) of DSA (Ohira et 
al. 2010a; 2010b). We can safely measure £^max(0 and r](t) 
only in such SNRs, since the effect of the secondary shocks 
on the spectral slope would not be dominant there. 
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